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TIME  DOMAIN  FREQUENCY  STABILITY 
CALCULATED  FROM  THE  FREQUENCY  DOMAIN 

DESCRIPTION; 

Use  of  the  SIGINT  Software  Package  to 
Calculate  Time  Domain  Frequency  Stability  From 
the  Frequency  Domain 

F.  L.  Walls  John  Gary  Abbie  O’ Gallagher 
Roland  Sweet  Linda  Sweet 

National  Institute  of  Standards  and  Technology 
Boulder,  Colorado  80303-3328 


We  describe  the  use  of  SIGINT,  an  interactive  software  package  devel- 
oped by  the  National  Institute  of  Standards  and  Technology,  which  facil- 
itates the  calculation  of  time  domain  frequency  stability  from  frequency 
domain  data  as  a function  of  measuring  time  in  terms  of  either  the  Al- 
lan variance,  o-y{T)]  the  modified  AUan  variance,  mod(7^(T);  or  (7x(r)  = 
(r/>/3)niodcry(T).  Except  for  the  graphic  output,  the  code  is  written  in 
standard  FORTRAN  77  and  runs  on  AT  compatible  computers  that  have 
a math  co-processor.  It  also  runs  on  many  other  systems;  however,  calls  to 
an  available  graphics  library  will  need  to  be  substituted  for  those  that  are 
included  in  this  version.  The  program  uses  either  a user  defined  function 
for  the  input  noise  or  default  functions  that  describe  the  noise  types  com- 
monly found  in  oscillators,  amplifiers,  frequency  multiphers,  frequency 
dividers,  and  general  signal  processing  equipment  including  up  to  four 
coherent  bright  lines  in  the  noise  spectra.  These  default  functions  make 
it  simple  to  analyze  the  time  domain  frequency  stability  as  a function  of 
measuring  bandwidth  using  realistic  first-,  second-,  or  third-order  low-pass 
filters  or  the  simplified  infinitely  sharp  cutoff  parameter  fh-  The  default 
functions  are  also  set  up  to  examine  the  effect  of  various  servo  parameters 
on  the  performance  of  a frequency  source  locked  to  a frequency  reference. 

Key  words:  Allan  variance;  frequency  lock  loop  analysis;  modified  Allan 
variance;  phase  lock  loop  analysis;  spectral  density  of  frequency  stability; 
time  domain  frequency  stability. 
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1 Introduction 


SIGINT  is  an  interactive  software  package  designed  to  facilitate  the  calculation  of  time 
domain  frequency  stability  from  data  in  the  frequency  domain.  The  Allan  variance 
(two-sample  variance),  cry[r)]  the  modified  Allan  variance,  mod(7y(r);  or  cTx(t)  = 
(t/\/3)  modcry(r)  can  be  calculated  as  functions  of  measurement  time.  Except  for 
the  graphic  output,  the  code  is  written  in  standard  FORTRAN  77  to  run  on  AT  and 
compatibles  that  have  a math  co-processor  (see  §4  for  details  on  portability).  The 
actual  output  of  the  software  is  the  square  root  of  the  variance.  The  equations  are 
[1,2,  3,4] 


cTy(nTo)  = 


(irfriToy 


rjh 

/ Sy{f) 

Jo 


df 


1 

2 


modcry(riTo)  = 


TI'^TT^Tq  1 JO 


+ 2 / XI cos(27r/fcro)  sin"^(7r /nro)  df 

k=i  ' 


sin'‘(7r/TiTo)  df 
Sy{f) 


P 


, and 


o-x(nro) 


(^)modCT„(nTo) 


where  the  measurement  time  is  given  by  nro  with  Tq  being  the  minimum  measurement 
time,  and  f is  the  Fourier  frequency  offset  from  the  carrier. 


The  input  parameters  are  expressed  in  terms  of  the  spectral  density  of  frequency 
fluctuations,  Sy{f).  Conversion  from  phase  noise  to  Sy{f)  is  simply 

Sy{f)  = 45*(/).and 

Sjf)  = ^Sy{f), 

where  uq  is  the  carrier  frequency. 

You  may  choose  to  supply  your  own  function  5y(/),  or  to  use  the  functions  which 
are  built  into  the  package.  Use  the  parameter  SELSY  (see  below)  to  communicate 
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this  choice  to  the  package.  If  you  use  your  own  function,  you  must  modify  the  empty 
FORTRAN  function,  SYF(f),  within  the  package.  The  built  in  function  has  the  form: 

‘^y(/)  = K{f)M{f)  ^ ^2/"^  + C3  + C4/  + Csf^ 

+ -1-  CsS(fms)  + CgS^fmg)}  H-  Cio  . 

The  values  of  the  Cx  are  stored  in  array  C.  Parameters  C\  through  C5  specify  the 
level  of  the  five  most  commonly  encountered  random  noise  types  found  in  oscillators, 
amplifiers,  frequency  multipliers,  frequency  synthesizers  and  general  signal  processing 
equipment.  Parameters  C&  through  Cg  specify  the  level  of  (up  to  four)  coherent  bright 
lines  in  the  frequency  noise  spectra.  That  is  for  i = 6,  7,  8,  9,  = Yrmsi-  Here 

S{fm)  = 1,  for  / = /m,  and 

S{fT7i)  = 0,  for/7^/m. 

Note  that  Yrmsi  ^Ho  be  expressed  as  where  (j)RMSi  is  the  RMS  value  of 

the  phase  modulation  of  Fourier  frequency  /m^. 

Servo  analysis,  such  as  the  examination  of  the  effect  of  various  servo  parameters  on 
the  locking  of  an  oscillator  to  a frequency  reference,  is  facilitated  by  choosing  the 
appropriate  form  of  K(f).  The  noise  type  assumed  for  the  reference  source  via  Cio,  is 
white  frequency  modulation.  An  example  of  this  is  given  in  the  Appendix.  In  some 
special  servo  cases,  it  may  be  necessary  to  divide  Sy{f)  into  two  segments  and  run 
them  as  separate  cases  , or  you  can  enter  your  own  form  for  5y(/). 

The  filter  function  K{f)  must  be  chosen  from  the  following  four  functions  by  setting 
the  value  of  SELK  (SELK=0,  1,  2 and  3 respectively). 

K{f)  = 1, 


The  coefficients  K{  must  be  input  in  the  array  CK. 
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The  upper  cutoff  frequency  of  the  integral  can  be  treated  as  infinitely  sharp  by  setting 
fh.  equal  to  the  equivalent  noise  bandwidth  of  the  time  domain  configuration  simulated 
by  the  calculations,  or  it  can  be  treated  more  realistically  using  the  appropriate 
first-,  second-,  or  third-order  low-pass  filter  function.  This  option  is  implemented  by 
choosing  the  appropriate  constants  for  Mi,  M2,  and  M3,  and  setting  fh  to  a sufficiently 
large  value  that  the  termination  of  the  integration  at  this  value  does  not  significantly 
bias  the  results. 


M(/)  =:  (l+Ml/(l+M2/)(l+M3/))^ 

where  any  or  all  of  the  M{  can  be  0.  These  Mi  are  stored  in  the  CM  array. 

This  accommodates  the  most  commonly  encountered  types  of  random  noise  found 
in  oscillators,  amplifiers,  frequency  multipliers,  frequency  synthesizers  and  general 
signal  processing  equipment.  The  rest  of  the  parameters  are  described  in  the  Input 
Parameter  section,  below. 


2 Running  the  Program 

The  SIGINT  package  of  routines  is  designed  for  interactive  use.  The  user  controls 
the  code  by  answering  questions  which  are  displayed  on  the  screen.  The  first  of  these 
questions  asks  whether  to  run  a new  case  or  read  output  from  a file.  (Output  to  be 
read  must  have  been  written  by  SIGINT  so  if  this  is  the  first  time  you  have  run  the 
code  there  will  be  no  output  available  to  read.)  If  you  choose  to  read  output  from  a 
file  the  code  will  ask  for  the  file  name,  read  in  that  entire  file  and  then  display  the 
“output  menu”  (see  below).  If  you  choose  to  run  a new  case  the  package  will  display 
the  values  of  the  parameters  that  define  the  default  case.  It  then  presents  a menu 
which  asks  whether  you  want  to  run  that  case,  be  prompted  to  input  an  entire  new 
set  of  parameters,  or  modify  directly  the  current  values  of  some  parameters.  This  will 
be  referred  to  later  as  the  “input  menu”. 
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2.1  Input  Modes 


In  either  of  the  modes  described  below,  values  may  be  entered  in  either  fixed-  or 
floating-point  format.  If  floating-point  is  used,  the  decimal  point  is  necessary;  for 
example,  you  can  enter  2.e-3  but  not  2e-3.  In  some  cases,  when  an  integer  is  being 
entered,  the  decimal  point  must  be  omitted.  In  all  cases,  the  program  recovers  from 
incorrect  input,  outputs  a message  attempting  to  diagnose  the  problem,  and  allows 
re-entry  of  the  value. 


2.1.1  The  “prompting”  input  mode 

If  you  elect  to  be  prompted  for  the  values  of  the  variables,  the  program  will  lead  you 
through  a series  of  questions  and  answers  by  which  a new  case  will  be  defined.  An 
advantage  of  this  mode  is  that  it  is  not  necessary  to  know  the  names  of  the  variables 
which  the  program  uses. 


2.1.2  The  “direct”  input  mode 

If,  instead,  the  “direct”  mode  of  input  is  chosen,  values  are  changed  by  typing  one 
or  more  lines  containing  variable  names  followed  by  “=”  and  the  new  values  of  the 
variable.  Different  entries  can  be  separated  by  a comma  or  blanks.  The  last  line  must 
be  terminated  by  a “$”  or  a For  example,  you  might  type  the  following  lines: 

NRANGE=2,  NLOW  = 10 
SELSY=1  CK=1. ,2.,3.  $ 

This  mode  will  usually  be  faster  than  the  “prompting”  mode. 

Notice  that  if  you  want  to  get  out  of  this  mode  without  changing  anything,  you  can 
simply  type  a “$”  or  a “;”  . 

When  you  finish  with  either  mode  of  input,  the  current  set  of  parameters  will  be 
displayed  and  you  will  be  returned  to  the  input  menu.  At  this  point  you  can  either 
run  the  case  or  further  modify  the  parameters  by  again  entering  one  of  the  input 
modes.  You  can  continue  to  revise  and  examine  the  parameters  until  you  are  satisfied 
with  them  and  then  run  the  case. 
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2.2  Output  Modes 


When  the  integration  is  completed  or  output  has  been  read  from  a file,  you  will  be 
presented  with  another  menu  which  will  be  called  the  “output  menu”.  This  menu  asks 
whether  to  print  results,  plot  them,  compute  another  integral,  read  output  from  a file 
or  quit.  As  the  program  computes  an  integral,  it  saves  the  defining  parameters  and 
the  results  in  internal  data  structures.  The  results  are  in  the  form  of  a table  of  values 
oi  T (t  = TiTo)  and  corresponding  values  of  the  square  root  of  the  integral  (variance). 
Any  results  you  read  from  output  files  will  be  in  the  same  form.  So  each  time  you 
come  to  the  output  menu,  the  results  of  all  of  the  cases  you  have  done  (or  read) 
during  the  run  are  available  to  be  printed  and/or  plotted.  If  you  choose  any  of  the 
output  modes,  and  if  any  of  the  cases  you  choose  to  output  are  modcTy  cases,  you  will 
be  given  the  option  of  having  those  values  converted  to  values  as  they  are  output. 
The  values  in  the  internal  arrays  will  remain  unchanged. 

2.2.1  Printing 

On  AT  compatible  computers  you  can  either  print  on  the  screen  or  on  your  printer 
(see  Portability  section  for  behavior  on  other  systems).  Note  that  if  you  are  working 
on  a system  that  has  no  printer  attached  and  you  choose  the  printer  option,  it  is 
likely  that  the  system  will  hang,  up  and  you  may  lose  your  results. 


2.2.2  Plotting 

Plotting  will  work  only  on  AT  compatible  computers.  On  each  frame  it  is  possible  to 
plot  either  one  or  several  cases.  A hard  copy  of  the  plot  can  be  obtained  if  you  have 
a memory- resident  utility  that  allows  you  to  use  the  print  screen  key  to  transfer  the 
graphics  screen  image  to  your  printer.  The  plot  package  will  plot  directly  to  certain 
printers.  Examples  of  the  plots  are  found  in  figures  1-9.  Additional  details  can  be 
obtained  by  contacting  the  Apphed  and  Computational  Mathematics  Division  of  the 
Computing  and  Apphed  Mathematics  Laboratory,  NIST,  325  Broadway,  Boulder, 
Colorado  80303-3328.  Attention:  Abbie  O’GaUagher. 
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2.2.3  Quitting  and  Saving  Results 

When  you  choose  to  quit,  you  will  be  asked  if  you  want  to  save  the  results.  If  you  say 
yes,  you  will  be  asked  for  the  name  of  a file  to  contain  these  results.  If  a file  exists 
under  that  name,  you  will  be  given  a chance  to  choose  another  file  name  or  write  over 
the  file.  If  there  are  moday  cases  present,  you  will  be  asked  whether  to  output  them 
as  (Ji  cases.  (Note  that  if  you  write  cases  out  as  you  will  not  be  able  to  recover 
the  modcTy  values  later.)  Finally  the  defining  parameters,  along  with  the  results  of 
each  of  the  cases  you  have  done  or  read  during  that  run,  will  be  saved  in  a disk  file 
on  the  default  disk  and  directory  under  the  name  you  have  specified.  The  next  time 
you  run  SIGINT  you  can  read  in  this  file  and  display  those  results. 


2.3  Input  Parameters 

The  following  is  a list  of  the  input  variables  that  must  be  set  in  order  to  compute 
an  integral.  The  integral  can  be  computed  for  a single  value  of  r (r  = nro),  or  for  a 
sequence  of  values  of  r.  In  the  latter  case,  the  integral  can  be  plotted  as  a function 
of  T. 


INTGRL:  Integer  which  selects  the  integral  to  be  computed.  INTGRL  = 1 for  the 

Gy  integral;  INTGRL  = 2 for  the  modo-y  integral;  INTGRL  = 3 for  the 
Gx  integral.  The  default  value  is  INTGRL  = 1. 


NRANGE:  Integer  which  determines  the  range  of  values  of  the  parameter  n for 

which  the  integral  is  computed.  Set  NRANGE  = 1 to  compute  the 
integral  for  a single  value  of  n (n  = NLOW).  Set  NRANGE  = 2 to 
compute  the  integral  for  the  sequence  of  values  obtained  by  starting 
with  NLOW  and  doubling  the  current  value  to  obtain  the  next  one  until 
reaching  the  last  such  value  not  exceeding  NHIGH.  Thus  if  NRANGE 
= 2,  NLOW  = 1,  and  NHIGH  = 20,  then  the  integral  will  be  computed 
for  n=(l,  2,  4,  8,  16).  For  NRANGE  = 3,  the  integral  will  be  evaluated 
by  decades  with  five  points  per  decade.  Thus  if  NRANGE  = 3,  NLOW 
= 1,  NHIGH  = 100,  then  the  evaluation  will  be  done  for  n = (1,  2,  3, 
5,  7,  10,  20,  30,  50,  70,  100).  The  default  value  is  NRANGE  = 3. 

NLOW:  The  lower  limit  of  the  range  of  values  of  n for  which  the  integral  is 
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NHIGH: 


FH: 

SELSY: 


SELK: 


C: 


FM: 


CK: 


CM: 


TAUO: 


computed  (see  the  discussion  of  NRANGE).  The  default  value  is  NLOW 

= 1. 

The  upper  limit  of  the  range  of  values  of  n for  which  the  integral  is 
computed.  The  default  value  is  NHIGH  = 1000. 

The  upper  limit,  f’^,  of  the  integral.  The  default  value  is  FH  = 3. 

Integer  to  select  the  function  Sy{f).  If  SELSY  = 1,  you  must  supply 
a double-precision  function  named  SYF(f),  with  double-precision  argu- 
ment f.  (You  must  replace  the  empty  FORTRAN  function  SYF(f)  that 
is  provided.)  If  SELSY  = 2,  then  you  must  select  one  of  four  built-in 
functions  by  giving  values  of  the  input  parameter  SELK,  and  the  arrays 
C,  CK,  and  CM  as  given  in  the  discussion  of  Sy{f)  above.  The  default 
value  of  SELSY  is  SELSY  = 2. 

Integer  parameter  which  selects  the  function  K(f).  See  discussion  above. 
The  default  value  is  SELK  = 0. 

A double-precision  array  of  with  10  elements  containing  the  coefficients 
of  the  function  Sy{f).  The  default  values  are  C — 2.e-24,  0,  0,  0,  0,  0, 
0,  0,  0,  0. 

A double-precision  array  with  4 elements,  FMs  through  FMg,  specify- 
ing the  modulation  frequencies,  /m*,  corresponding  to  the  sixth,  sev- 
enth, eighth  and  ninth  elements  of  the  C array.  The  default  values  are 
0.,  0.,  0.,  0.  However  note  that  if  Ci  is  nonzero,  FM{  must  be  greater 
than  zero.  Also,  for  Ci  to  contribute  to  results,  FMi  must  be  less  than 
FH. 

A double-precision  array  with  3 elements  containing  the  coefficients  of 
the  function  K(f).  The  default  values  are  CK  = 10,  40,  100.  Note, 
however,  that  these  values  do  not  figure  in  the  calculation  when  SELK 
has  it’s  default  value  of  0. 

A double-precision  array  with  3 elements  containing  the  coefficients  of 
the  function  M(f).  The  default  values  are  CM  = 0 ,0,  0. 

A double-precision  scalar  parameter  of  the  integral.  The  default  value 
is  TAUO  = 1. 
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2.4  Internal  Constants 


It  is  not  necessary  to  be  aware  of  the  constants  described  in  this  section  unless  you 
want  to  alter  the  behavior  of  the  code.  In  that  case  it  may  be  necessary  to  change  the 
values  to  which  these  constants  are  set.  This  may  be  done  by  editing  the  FORTRAN 
to  change  the  assignment  statements  that  set  these  values.  Recompiling  will  also  be 
necessary.  The  pertinent  assignment  statements  can  be  found  at  the  beginning  of  the 
main  program. 


EPSREL:  The  relative  error  tolerance  for  QUADPACK  routines.  It  is  set  at 

2*  10-^ 

EPSTRN:  The  truncation  error  tolerance  used  by  the  MDSIGY  routine  to  deter- 

mine, for  a given  value  of  the  variable  of  integration,  whether  the  esti- 
mate of  the  integral  over  the  remainder  of  the  interval  is  small  enough 
to  indicate  that  the  remainder  can  be  ignored.  See  discussion  of  the 
numerical  method,  below.  EPSTRN  is  set  at  2 * 10“^. 


NCYCLE:  The  number  of  cycles  of  the  numerator  of  Sy{f)  on  each  side  of  a sin- 

gularity which  are  evaluated  by  the  non-oscillatory  quadrature  routine 
(see  section  3).  NCYCLE  is  set  to  7.1. 


IDB:  Switch  to  control  debug  printout  from  the  routines  SIGMAY  and  MD- 

SIGY. It  is  set  to  0.  Set  it  to  1 to  get  debug  printing. 


3 The  Numerical  Method 


The  sigma-y  integral, 

1. 
3 


can  be  treated  as  an  oscillatory  integral  except  in  the  neighborhood  of  the  origin. 
Therefore,  the  integral  is  computed  by  using  one  approximation  in  the  neighborhood 
of  the  origin  and  a second  method  away  from  the  origin.  The  bright  line  terms. 


(7y{nTo)  = 


L 
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those  which  involve  the  coefficients  Ce>  through  Cg,  can  be  calculated  directly.  In  the 
following  discussion,  then,  let 


Then 


where 


and 


C\f  ^ + C2/  ^ + C3  H-  C4/  + CsP 


— [2(A  + -^2  + ^3)]^ ) 

h = rF,{f)df, 

Jo 

I2  = [ F2{f)  sm^{TT  fnro)df, 

J a 

^ 2Ci  sin^^TT  fruinTo) 


Fiif) 

m) 


= Syif) 


sin^(7r/nro) 

(7r/nro)2 


Sy{f) 

(7r/nro)2 ' 


M(f)  and  K(f)  are  as  described  in  the  introduction,  a is  chosen  so  that  the  first 
integral  is  taken  over  NCYCLE/2  cycles  of  the  sine  function,  provided  r = nro 
is  large  enough.  Otherwise  it  is  taken  over  2/3  of  the  interval  fh.  That  is,  a = 
mzmmum(2//i/3,  NCYCLE/r).  If  the  function  Fi  grows  more  slowly  than  at 
/ = 0,  then  the  first  integral  exists.  We  assume  that  Sy{f)  grows  no  more  rapidly 
than  and  thus  the  integrand  in  this  first  integral  has  a removable  singularity  and 
can  be  computed  using  a general  purpose  quadrature  routine  from  the  subroutine 
library.  We  use  the  DQAGE  routine  from  the  QUADPACK  [7]  hbrary. 


The  second  integral  is  oscillatory.  It  can  be  transformed  into  a canonical  form  by  use 
of  the  trigonometric  substitution, 


3 1 1 

sin^^TTfr)  = cos(27r/r)  + - cos(47t/t), 

8 2 8 
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and  can  therefore  be  written  as  the  sum  of  three  integrals, 


I2 


^21  + -^22 


'23- 


Since  a is  chosen  large  enough  to  avoid  the  singularity  at  the  origin,  the  first  of 
these  three  integrals  can  be  computed  with  the  standard  quadrature  routine  DQAGE. 
The  remaining  two  integrals  are  oscillatory  and  are  therefore  computed  using  the 
DQAWOE  routine,  also  from  QUADPACK.  DQAWOE  is  designed  for  integrals  of 
the  form  / f(x)sin{ujx)dx  or  / f{x)cos{ux)dx  , where  uj  can  be  very  large.  These 
sigma  integrals  may  have  values  of  nro  which  easily  exceed  10^  (here  titq  corresponds 
to  Lj). 

In  all  cases  we  ask  that  the  QUADPACK  routines  compute  the  integrals  with  a 
relative  error  less  than  2 * 10“^.  This  may  be  wasteful,  since  one  of  these  integrals 
may  dominate  the  remainder,  however  we  could  not  be  certain  of  the  dominance  in 
all  cases. 


Extensive  tests  of  the  numerical  calculation  for  o‘y(r)  have  been  performed  for  all  the 
common  noise  types,  C\  — C5  [1,  2,  3].  The  results  agree  with  the  asymptotic  form 
(27r/fcT  ^ 1)  to  better  than  11%  as  long  as  2tt fhT  is  larger  than  6.3  and  0.1%  for 
2TrfhT  > 62.8.  We  believe  that  the  numerical  technique  yields  more  accurate  results 
since  the  analytical  calculations  are  valid  only  for  2'KfhT  ^ 1.  The  asymptotic  forms 
are  considered  in  more  detail  in  Figures  1-7  and  Tables  I and  II. 


The  mod  sigma-y  integral, 


modo-y(nro)  = 4 ^2  2 | ^ / sin^(7r/nro)  d/ 

71  TT^Tq  I JQ  J 

+ 2 [ y'(n  - cos(27t/A;to)  sin'^(7r/nTo)  df 


is  much  more  difficult  to  approximate  because  the  sum  over  k includes  many  terms 
when  n is  large.  The  sum  over  k is  eliminated  by  using  the  following  identity. 


n— 1 

y^(7i  — k)  COS  6k 
k=l 


n 1 sin^  ~ 

2 2 sin^  I 


This  identity  is  derived  by  summing  a geometric  series  composed  of  complex  expo- 
nentials. Since  the  real  part  of  a complex  exponential  is  the  cosine  term,  this  series 
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reduces  to  the  expression  on  the  left,  above.  This  can  be  summed  explicitly,  yielding 
the  final  formula.  Application  of  this  identity  to  the  original  integral  leads  to  the 
following  form  of  the  integral; 


modcTy(nTo) 


2 /-/h  5y(/)sin^(7rron/) 

n'^TT^TQ  7o  sin^(7TTo/) 


With  this  simplification  modcry(nro)  can  be  broken  up  and  evaluated  in  a way  that 
is  analogous  to  that  described  above  for  the  dyinro)  integral.  However  modcry(nro) 
presents  the  additional  problem  that  it  has  a singularity  for  each  integral  value  of 
To/  instead  of  just  one  singularity  at  / = 0.  It  is  therefore  necessary  to  break  up  the 
interval  of  integration  into  fh  * Tq  subintervals  and  sum  the  values  of  the  integral  over 
each  of  these  smaller  intervals.  If  fh  is  large,  this  is  a large  number  of  subintervals.  But 
in  many  cases,  the  later  subintervals  contribute  negligibly  to  the  value  of  the  integral. 
For  this  reason  an  estimate  of  the  error  produced  by  neglecting  these  contributions 
is  made  every  few  subintervals.  If  the  ratio  of  this  error  to  the  computed  value  of  the 
integral  up  to  that  point  is  small  (less  than  EPSTRN),  then  no  further  calculation  is 
done  and  the  value  returned  by  the  program  as  modo-y(nTo)  is  the  square  root  of  the 
sum  of  the  integrations  done  over  the  subintervals  up  to  that  point. 

As  in  the  cjy  case,  the  bright  line  terms  (those  which  involve  the  coefficients  Ce  through 
Cg)  are  calculated  directly.  So,  in  fact,  the  integral  becomes: 


moday(riTo) 


2 Sy{f)sm%TrTonf) 

Jo  P sin^(7rro/) 

^ 2 (7t  sin®  ( 7T  frui  titq  ) 

^ (tt  frrii  nToYri^ sin^{'K  frrii  Tq)M{  frrii  )K{  fm^ ) 

Cif-^  + C^f-^  + C3  + C,f  -f  CsP 
K{f)M{f)  ^ 


3 

j 


with  M(f)  and  K(f)  as  described  in  the  introduction. 


Extensive  tests  of  the  numerical  calculation  for  modcTy(T)  have  been  carried  out  for 
the  common  types  of  noise  {Ci  — Cs)  and  compared  to  the  results  of  [3],  obtained 
analytically  and  the  results  of  [4]  obtained  by  computer  simulated  noise.  This  has  been 
done  by  comparing  R(n)  the  ratio  of  modcTy(r)  to  as  a function  of  the  number  of 
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samples  averaged  together  for  modcTy(T),  obtained  by  the  three  methods.  The  results 
obtained  using  SIGINT  for  R(n)  are  shown  in  figure  1 and  table  1.  This  agrees  exactly 
with  the  results  of  both  [3]  and  [4]  for  a = —2,  a = 0 and  a = 2 (ignoring  the  obvious 
typographical  errors).  The  results  for  a = 1 agree  for  27t//iTo  = — the  only  case 

addressed  by  [3].  Also  shown  are  our  results  for  a = 1 and  2t^ = 3,  10,  100.  Our 
results  differ  considerably  from  [3]  for  a = — 1,  but  agree  with  the  asymptotic  value 
of  [4].  Since  we  have  obtained  excellent  agreement  for  so  many  other  cases  using  the 
same  numerical  code,  we  believe  that  the  results  presented  in  figure  1 and  table  1 are 
more  accurate  than  the  earlier  work. 


4 Portability 

The  SIGINT  code  is  standard  FORTRAN  77  and  is  portable  with  the  following  ex- 
ceptions: 

• The  graphics  library  used  is  David  Kahaner’s^  GRAPH. LIB,  which  runs  only 
on  AT-compatible  computers  with  one  of  the  following  display  adapters: 

— AT-compatible  Color  Graphics  Adapter^ 

— AT-compatible  Enhanced  Graphics  Adapter 

- AT-compatible  Professional  Graphics  Adapter  in  CGA  emulation  mode 
— AT-compatible  PS/2  Video  Gate  Array  in  EGA  emulation  mode 
— Hercules 
— Video- 7 

When  moving  the  code  to  another  system,  it  will  be  necessary  either  to  “com- 
ment out”  the  calls  to  the  plot  routines  and  run  the  code  without  graphic  output 
or  to  change  those  graphics  calls,  all  of  which  are  in  SUBROUTINE  PLOT,  to 
caJl  an  available  graphics  library. 

^Scientific  Computing  Division,  National  Institute  of  Standards  and  Technology,  Gaithersburg, 
MD  20899 

^Identification  of  some  commercial  materials  has  been  necessary  in  this  report.  In  no  case  does 
such  identification  imply  recommendation  or  endorsement  by  the  National  Institute  of  Standards  and 
Technology,  nor  does  it  imply  that  the  materials  are  necessarily  the  best  available  for  the  purpose. 
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• On  any  system,  numerical  results  can  be  printed  on  the  screen  but  if  the  option 
of  using  the  printer  is  chosen  on  a system  other  than  an  AT  compatible,  the 
results  will  instead  be  written  to  a file  called  “PRN”.  (“PRN”  is  the  name  DOS 
uses  for  the  printer.)  Results  will  not  accumulate  in  this  file.  Only  the  last 
printing  will  be  there  when  the  program  terminates. 

• The  code  is  written  in  double  precision,  which  may  be  unnecessary  on  machines 
with  larger  word  sizes. 

• On  UNIX  systems,  carriage  control  characters  may  appear  on  the  screen  instead 
of  affecting  the  output  as  intended. 

• The  VAX  reports  underflow  conditions  but  terminates  normally  and  the  results 
are  unaffected. 


The  entire  code  has  been  run  extensively  on  AT-compatible  computers  using  RM/ FOR- 
TRAN. In  addition,  except  for  the  graphic  output,  it  has  been  tested  on  a VAX 
11/785,  SUN  3/180,  MASSCOMP  MC5500-PEP,  and  Cyber  840  and  855.  It  is  ex- 
pected that  it  will  run  on  most  other  systems.  If  problems  are  encountered  in  moving 
it  to  another  system,  an  attempt  will  be  made,  subject  to  the  availability  of  resources, 
to  overcome  them.  For  such  assistance,  please  contact  the  Applied  and  Computa- 
tional Mathematics  Division,  NIST,  325  Broadway,  Boulder,  Colorado  80303-3328. 
Attention:  Abbie  O’Gallagher. 
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Figure  1:  CTy{T)  versus  r for  the  five  common  power-law  noise  types  in  the 
limit  that  2t^  fhT  is  large  compared  to  1 and  an  infinitely  sharp  filter  is  used. 
Curve  a is  for  random- walk  frequency  modulation,  Sy{f)  = Curve 

b is  for  flicker  frequency  modulation,  Sy{f)  = Curve  c is  for  white 

frequency  modulation,  Sy{f)  = ho.  Curve  d is  for  flicker  phase  modulation, 
Sy{f)  = hif  and  fh  = 16  Hz.  Curve  e is  for  white  phase  modulation, 
Sy{f)  = h2p  and  fh  — 16  Hz.  /12  = /ii  = /lo  = ^-1  = ^-2  = 2x10“^'* 
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Figure  2:  a-y(r)  for  white  phase  modulation  (a  = 2 and  /12  = 2il0~^^)  as  a 
function  of  measurement  time,  r,  and  measurement  bandwidth,  f^.  Curves 
a,  b,  and  c have  an  infinitely  sharp  filter  with  width,  fh  = 16  Hz,  fh  = 0.016 
Hz,  fh  = 0.0016  Hz  respectively.  Curves  d and  e have  a single  pole  filter 
width,  fh  = 0.016  Hz  and  fh  = 0.0016  respectively. 
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Figure  3:  cry(r)  for  flicker  phase  frequency  modulation  (a  = 1 and  hi  = 
2x10“^"^)  as  a function  of  measurement  time,  r,  and  measurement  bandwidth, 
fh.  Curves  a,  b and  c have  an  infinitely  sharp  Alter  with  width,  fh  = 16  Hz, 
= 0.016  Hz,  fh,  = 0.0016  Hz  respectively.  Curves  d and  e have  a single 
pole  Alter  width,  fh  = 0.016  Hz  and  fh  = 0.0016  respectively. 
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Figure  4:  CTy{r)  for  white  frequency  modulation  (a  = 0 and  ho  = 2x10“^'^)  as 
a function  of  measurement  time,  r,  and  measurement  bandwidth,  fh-  Curves 
a,  b and  c have  an  infinitely  sharp  filter  with  width,  fh  = 16  Hz,  fh  = 0.016 
Hz,  fh  = 0.0016  Hz  respectively.  Curves  d and  e have  a single  pole  filter 
width,  fh  = 0.016  Hz  and  fh  = 0.0016  respectively. 
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Figure  5:  cry(r)  for  flicker  frequency  modulation  (a  = —1  and  /i_i  = 2xl0~^'^) 
as  a function  of  measurement  time,  r,  and  measurement  bandwidth,  fh- 
Curves  a,  b and  c have  an  infinitely  sharp  filter  with  width,  fh  = 16  Hz, 
= 0.016  Hz,  fh  = 0.0016  Hz  respectively.  Curves  d and  e have  a single 
pole  filter  width,  fh  = 0.016  Hz  and  fh  = 0.0016  respectively. 
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Figure  6:  CTy(r)  for  random-walk  frequency  modulation  (a  = —2  and  h_2  = 
2ccl0~^'^)  as  a function  of  measurement  time  and  measurement  bandwidth, 
fh  for  7T.  Curves  a,  b and  c have  an  infinitely  sharp  filter  with  width,  //^  = 16 
Hz,  fh  = 0.016  Hz,  fh  = 0.0016  Hz  respectively.  Curves  d and  e have  a single 
pole  filter  width,  fh  = 0.016  Hz  and  fh  = 0.0016  respectively. 
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n,  Number  of  Samples  Averaged 


Figure  7:  Ratio  of  to  moday(tau)  as  a function  of  n,  the  number  of 

points  averaged  to  obtain  mo dcTy{tau).  The  measurement  time  r = nro, 
where  tq  is  the  minimum  data  interval. 
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Figure  8:  Curve  a shows  (7y(nro)  for  white  phase  modulation  (o  = 2,/i2  = 
2x10“^'*,  fh  = 16  Hz  and  coherent  frequency  modulation  Cq  = — 10~^^ 

at  a frequency  of  fme  = 6 Hz.  Curve  b shows  <7y(nro)  for  the  same  noise 
parameters  as  curve  a with  Ce  = 0.  The  integrals  were  evaluated  at  n = 1, 
2,  3,  5,  7,  10,  20,  ...10,000  and  tq  = 0.008  ms. 
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Figure  9:  Curve  a shows  mod  cTy(nTo)  for  white  phase  modulation  (a  = 
2, /i2  = 2x10“^"^,//,  = 16  Hz  and  coherent  frequency  modulation  Cq  = 
^RMS  — 3,t  a frequency  of  fm^  = 6 Hz.  Curve  b shows  mod  cTy(nro) 

for  the  same  noise  parameters  as  curve  a and  Cq  = 0.  The  integrals  were 
evaluated  at  n = 1,  2,  3,  5,  7,  10,  20,  ...10,000  and  tq  = 0.008  ms. 
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TABLE  1 Ratio  of  mod(7y(T)  to  <7y(r)  versus  n in  the  limit  that  > 1 , for 

common  power-law  noise  types  Sy(f)  = n is  the  number  of  time 

or  phase  samples  averaged  to  obtain  mod(7'^{T  = nro)  where  tq  is  the 
minimum  sample  time,  and  Uh  is  2tt  times  the  measurement  bandwidth 

Ih. 


R(n)  = 

modOy(T) 

VS. 

n t = 

^^0 

o^/t) 

n 

a = -2 

a = - 1 

1 

0 

II  j 
C j 

o 

II 

a 

Tq  = 10 

= +1 

t^»-o  = 100 

o 

II 

o 

f 

a = -►2 

1 

1.000 

1.000 

1.000 

1 .000 

1 .000 

1.000 

1.000 

1.000 

2 

0.859 

0.738 

0.616 

0.568 

0.543 

0.525 

0.504 

0.500 

3 

0.840 

0 701 

0.551 

0 481 

0.418 

0.384 

0.355 

0.330 

4 

0.831 

0.681 

0.530 

0.405 

0.359 

0.317 

0.284 

0 250 

5 

0.830 

0.684 

0.517 

0.386 

0.324 

0.279 

0.241 

0.200 

6 

0.828 

0.681 

0.514 

0.349 

0.301 

0.251 

0.214 

0.167 

7 

0 827 

0.679 

0.507 

0.343 

0.283 

0.235 

0.195 

0.143 

8 

0.827 

0 678 

0.506 

0.319 

0.271 

0.219 

0.180 

0.125 

10 

0.826 

0.677 

0.504 

0.299 

0.253 

0.203 

0.160 

0.100 

14 

0.826 

0.675 

0.502 

0.274 

0.230 

0.179 

0.137 

0.0714 

20 

0.825 

0.675 

0.501 

0.253 

0.210 

0.163 

0.119 

0.0500 

30 

0.825 

0.675 

0.500 

0.233 

0.194 

0.348 

0.106 

0.0333 

50 

0.825 

0.675 

0 500 

0.210 

0.176 

0.134 

0.0938 

0.0200 

100  0.825 

0.675 

0.500 

0.186 

0.159 

0.121 

0.0837 

0.0100 

Liaic  0.825 

0.675 

0.500 

f 

3 37  _ ] 

1/n 

1 1 

.04  -t-3  r J 
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TABLE  2 Asymptotic  forms  of  for  various  power-law  noise  types  and  two 

filter  types.  Note:  u’/,/27t  = f^^  is  the  measurement  system  bandwidth 
- often  caJled  the  high-frequency  cutoff.  In  = loge. 


Nane  of  Noise 

a 

s,(f) 

o*(r) 

u^r»l 

Infinite  Sharp 
Filter 

Whr»l 

Single  Pole 
Filter 

<^T«l 

Infinite  Sharp 
Filter 

u>hr«l 

Single  Pole 
Filter 

White  Phase 

2 

3ff,  h2 

K f ^ ■ ■ ■■■ 

3fhhz 

2/5n^  fh'^"h2 

fh'h2 

(2ir)2f2 

(2)r)2r2 

2r 

Flicker  Phase 

1 

(1.038  + 31n(whf))hi 

(31n(whf 

))hi 

2fh^(ln(2))hi 

hi  f 

(2ir)2f2 

(2>r)2 

r2 

White  Frequency 

0 

ho 

ho  

2r 

ho 

2r 

2/3>r^fh^  r^h(3 

2/3jT^fh^rho 

Flicker  Frequency 

-1 

h.if-^  2(ln(2))h_i 

2(ln(2))h 

-1  ir2f„2^2^i  ^ 

8ir^  fi,  ^ h.  1 

Random -Walk 
Frequency 

-2 

2w^  rh. 2 

u 2 

2jf2  rh-j 

2*2  fhr2h_2 

2»2f^r2h.2 

rk.  2 ^ 

3 

3 
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Stability  of  Frequency  Locked  Loops 
Fred  L.  Walls 

National  Institute  of  Standards  and  Technology 
Boulder,  Colorado  80303 


1 ■ Introduction 


Passive  frequency  standards  are  characterized  by  the  use  of  a reference 
resonance  to  stabilize  the  frequency  of  an  external  probe  oscillator.  A 
common  configuration  is  shown  in  Fig.  1 [1-6].  The  probe  oscillator 

generally  has  phase  modulation  imposed  on  the  carrier  in  order  to 
interrogate  the  resonance  with  a minimum  of  offset.  The  resulting 

amplitude  modulation  is  demodulated  to  yield  an  error  curve  that  is 
essentially  the  derivative  of  the  resonance.  Although  Fig.  1 shows  a 

transmission  system,  similar  schemes  are  sometimes  possible  in  reflection 
[5].  The  error  signal  from  the  demodulator  is  used  to  steer  the  probe 
signal  toward  the  center  of  the  resonance  line  [1-9].  For  analysis  times 
longer  than  one  period  of  the  modulation  cycle,  and  under  the  condition 
that  the  probe  oscillator  wanders  less  than  the  half  width  of  the  error 
curve  in  the  loop  attack  time,  we  can  treat  this  curve  as  approximately 
static[l-4].  Near  line  center  the  loop  error  voltage  Vp  at  the 
synchronous  detector,  is  approximately  = k (i/q®  - where  k is 

the  slope  of  the  error  curve,  i/r  is  the  resonance  frequency  of  the 
reference,  and  is  the  open  loop  frequency  of  the  probe  oscillator  and 

is  the  detector  noise.  If  we  now  close  the  loop  with  gain  G(f),  it 

can  be  shown  that  the  spectral  density  of  fractional  frequency 

fluctuations  Sy(f)  for  the  probe  source  becomes 


MODULATION 

REF 


PROBE 

SOURCE 


Figure  1 . Generalized  block  diagram  of  a probe  source  locked  to  a 
reference  resonance. 


Contribution  of  the  U.S.  Government;  not  subject  to  copyright. 


26 


+ s/(f) 


(1) 


Sy(f)  = 


G(f)  ' 

2 

lG(f)  + lJ 

V 

R 

+ 


S/(f)^ 
G"(f)  . 


where  Sy°(f)  is  the  open-loop  spectral  density  of  fractional  frequency 
fluctuations  of  the  probe  source,  Sy^(f)  is  that  of  the  reference,  and 
Sy^(f)  is  that  of  the  detector  and  interrogation  noise  referred  to  the 
demodulator  output.  Cutler  [10]  has  pointed  out  that  noise  in  the  local 
oscillator  at  the  2nd  harmonic  of  the  modulation  frequency,  which  is 
usually  ignored,  causes  a time  varying  frequency  offset  that  is 
undistinguishable  from  reference  noise.  This  sets  the  lower  limit  to  the 
interrogation  noise  and  often  sets  the  lower  limit  of  the  noise 
performance  of  the  local  oscillator  necessary  not  to  degrade  the  overall 
performance.  The  magnitude  of  the  2nd  harmonic  noise  modulation  in 
radians/s  is  estimated  in  appendix  B of  [1]  to  be  k3  = 2^0/^  Sy  (Q/tt)  ) ^ ^ ^ , 
where  Q/(2n)  is  the  modulation  frequency.  This  leads  to  an  interrogatic" 
noise  term  which  is  of  order  [1,4,10]  Sy^(f)  = 1/(167t)  Sy°  (fi/Tr)  . For 
large  values  of  G(f),  Sy(f)  of  the  probe  source  reflects  that  of  the 
reference  plus  the  added  noise  of  the  detection  system  [1-6]. 

The  primary  goal  of  this  paper  is  to  investigate  the  effect  of  various 
realistic  forms  of  G(f)  on  the  spectral  density  of  frequency  and 
fractional- frequency  stability.  It  will  be  shown  that  mod  cTy(r)  [11,12] 
is  better  suited  than  the  traditional  two-sample  or  Allan  Variance 
[8],  for  evaluating  the  locked  performance  when  the  frequency  stabi  lity 
of  the  reference  is  much  greater  than  that  of  the  local  oscillator 
[7,8,11,12] . 

2.  The  Effect  of  Different  Forms  of  Servo  Gain 

Figure  2 shows  the  effect  of  locking  a probe  source  with  Sy°(f)  = 2 x 
10'^®/f^  + 1 X 10"  ^^/f  + 2 X 10"  ^°f^  (which  roughly  corresponds  to  that 
of  a low-noise  5 MHz  quartz  oscillator)  to  a reference  resonance  with 
Sy^(f)  of  2 X 10"^°  and  gy^(f)  = 5.6  x 10"^®.  Sy^(f)  is  the  estimated 
interrogation  noise  for  a modulation  frequency  of  47  Hz  [1,4,10].  Curves 
A,  B,  and  C show  the  effect  of  using  a first-,  second-,  or  third-  order 
loop  each  having  a loop  bandwidth  of  approximately  0.1  Hz  or  an  attack 
time  of  1.6  s[l].  The  solid  Curves  A,  B,  and  C of  Fig.  3 show  Oyir) 
calculated  for  curves  A,  B,  and  C from  Fig.  2 with  a noise  bandwidth  of  3 
Hz.  The  improvement  in  stability  of  the  local  oscillator  due  to  the 
servo  scales  roughly  as  r/r^  where  is  the  attack  time  and  r is  the 
measurement  time.  Considerable  insight  into  the  effect  of  various  gain 
stages  on  the  frequency  stability  can  be  obtained  by  considering 
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J / (2nu^T)^ 


(2) 


Ai/ 


u, 


as  a measure  of  fractional  time  or  frequency  stability  in  the  region  of 
measurement  times  from  about  1 to  10^  s where  the  effects  of  low 
frequency  divergence  can  often  be  ignored  [7,8].  The  squared  phase 
deviation  of  the  zero  crossings  is  given  approximately  by [8] 


(Sy  (f))/f2  df 

1 /27t^  r 


(3) 


Figure  2.  Sy(f)  of  the 
probe  source  locked  a. 
reference  resonance.  For  A 
G(f)  = (l/(10f) , for  B G(f) 
= (l/(10f))(l  + l/(40f)), 
for  C G(f)  = l/(10f))  (1  + 
l/(40f))(l  + l/(160f)). 
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Figure  3.  The  solid 
curves  show  a^ir)  and  the 
dashed  curves  show  mod 
o’y(T')  versus  r for  curves 
A,  B,  and  C of  Fig.  2. 

The  a,  b,  and  c points 
show  estimations  of  mod 
ay(r)  using  Eqs . 2 and  3 
with  f^  = f^o • Also  shown 
is  (r)  and  (r ) . 


1 10^  10“  10® 
MEASUREMENT  TIME  (S) 


for  a sample  time  r and  noise  bandwidth  fj^  . For  a noise  spectrum  which 
varies  as  Sy(f)=  K'f"  where  n > 2,  the  integral  is  dominated  by  the  high 
frequency  bandwidth  even  for  very  long  measurement  times.  The 
fractional-frequency  (or  time)  stability  given  by  Eq.  2 decreases  as  t~ ^ 
just  as  does  cTy(r),  due  to  the  increase  in  measurement  time  and  not  to  a 
decrease  in  the  value  of  A<f>^  . The  integral  in  Eq.  3 can  be  integrated  by 
parts  for  the  various  segments  of  Sy(f)  . This  makes  it  easier  to 
evaluate  and  optimize  the  performance  of  the  overall  system  than  by  using 
a process  based  on  ay(r). 

The  contribution  of  the  high  frequency  noise  can  be  greatly  reduced  by 
phase  averaging  the  data  points [9].  The  data  at  measurement  time  r = nr^ 
(where  is  the  data  interval)  is  obtained  by  averaging  the  n adjacent 
phase  points.  This  is  equivalent  to  using  mod  ay(r)  to  analyze  the 
data[ 9 , 11 , 12  ] . The  standard  expression  for  mod  o^(t)  contains  an 
enormous  number  of  terms  and  is  quite  laborious  to  compute [ 11- 12 ] . It 
[13-15]  has  been  pointed  out  that  mod  (7y(r)  can  be  reduced  to 


mod  (7y  ( r ) 


2 


n-  rl 


0 


«/ 


Sy  (f)sin®  (TTT^nf) 

df. 

f^  sin^  (iTT^  f ) 


(4) 


which  is  much  more  manageable.  Nevertheless  it  still  requires  numerical 
calculations  to  determine  which  segment  of  the  phase  noise  dominates  the 
integral.  We  can  estimate  mod  <7y(r)  from  Eqs.  2 and  3 by  using  ft^  = 
fj^g/n,  where  f^^^  is  the  hardware  bandwidth  of  the  measurement  system  and 
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r = nr^  is  the  measurement  time.  The  integral  for  in  Eq . 3 is  now 
substantially  reduced  for  measurement  times  large  enough  that  fho/^  is 
less  than  the  bandwidth  of  the  servo  system.  This  is  in  contrast  to  the 
original  calculation  for  Eq . 3 where  the  integral  must  always  increase 
with  r.  This  integral  can  also  be  divided  into  parts  and  integrated 
analytically.  The  fractional  frequency  stability  computed  for  A,  B,  and 
C of  Fig.  2 using  mod  a^{r)  are  shown  as  the  dashed  curves  A'  , B'  , and  C' 
of  Fig.  3.  The  points  labeled  a,  b,  c show  the  results  of  estimating  mod 
o^{t)  using  Eqs . 2 and  3 with  fj^  = f^^/n.  The  agreement  with  the  dashed 
curves  is  very  good.  The  calculations  for  and  mod  <7y(r)  are 
virtually  independent  of  using  an  upper  cutoff  frequency  for  the 
integration  or  a simple  low  pass  filter  of  the  same  bandwidth.  The  use 
of  mod  cTy(r)  reduces  by  a factor  of  about  100  the  time  necessary  to  reach 
the  performance  of  the  reference  plus  the  interrogation  noise,  which  in 
this  case  limits  the  performance  for  times  longer  than  100  s, 

3,  Discussion 


The  primary  utility  of  these  results  is  the  insight  into  the  origin  of 
the  major  contributions  to  a^{r)  and  mod  a^ir)  as  a function  of  the  servo 
gain  and  the  analysis  bandwidth,  and  the  effect  of  narrowing  the 
bandwidth  with  longer  measurement  times.  The  specific  examples 
illustrate  that  it  is  generally  necessary  to  use  a second-order  loop  to 
lock  the  probe  source  to  the  reference  resonance,  but  that  a third-order 
loop  offers  little  additional  improvement  when  drift  in  the  probe  is  not 
serious.  In  addition,  if  the  reference  resonance  is  substantially  more 
stable  than  the  probe  source,  it  can  be  very  useful  to  use  mod  a^(r)  to 
analyze  the  output  of  the  probe  source.  We  have  introduced  a simple 
measure  of  frequency  stability  that  is  easy  to  use  for  optimizing  the 
servo  and  the  analysis  system.  With  the  techniques  introduced  here  it  is 
possible  to  realize  fractional  frequency  stabilities  and  time  prediction 
of  the  local  oscillator  which  are  characteristic  of  the  reference 
resonance  in  a way  that  is  orders  of  magnitude  faster  than  those  using 
more  traditional  approaches. 
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